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1.  Introduction 

Modeling  time-series  with  linear  pole-zero  AutoRegressive-Moving  Average  (ARMA) 
models  has  numerous  applications  in  signal  processing.  This  problem  is  in  general  non  linear  and 
most  ARMA  modeling  techniques  are  iterative  in  nature.  The  Iterative  Prefiltering  (IP)  method 
has  the  advantage  of  computing  potential  non-minimum  phase  representations  which  may  be  useful 
in  time-domain  modeling.  The  original  IP  minimization  procedure  is  an  ill-conditioned  problem 
which  has  classically  been  solved  using  a  least-squares  approach.  This  work  presents  a  modification 
of  the  classical  IP  technique  in  which  the  least-squares  iteration  step  is  replaced  by  a  Total  Least 
Squares  (TLS)  step  to  take  advantage  of  the  statistical  properties  of  the  TLS  method.  Results  show 
that  improvements  in  the  modeling  performances  may  be  obtained  by  using  the  TLS-based  IP 
method  when  modeling  signals  distorted  by  white  Gaussian  noise. 

Section  2  presents  a  review  of  the  TLS  method,  and  illustrates  its  main  advantages  as 
compared  to  the  LS  method.  Next,  Section  3  presents  the  classical  Iterative  Prefiltering  method  and 
the  proposed  TLS-based  IP  method.   Finally,  conclusions  are  given  in  Section  4. 

2.  Total  Least  Squares  Problem 

2.1  Introduction 

The  Total  Least  Squares  method  (TLS)  has  been  introduced  recently  as  a  substitute  to  the 
classical  Least  Squares  (LS)  method  for  solving  overdetermined  linear  systems  of  equations  when 
all  data  involved  in  the  computation  are  corrupted  with  noise  (errors)  [1,2].  Consider  the  problem 
of  solving  the  overdetermined  system  of  equations  Ax  =  b,  where  A  has  dimension  m*n  and  b  has 
dimension  m*l.  The  LS  solution  xls  is  obtained  by  minimizing  the  Euclidean  norm  ||Ax-b||2. 
Solving  for  the  LS  solution  is  equivalent  to  solving 
the  following  minimization  problem: 

minimize  ^ejr    \\b-h!_\\2 
subject  to    bl  e  R(A) 


The  above  equation  is  satisfied  when  b'  is  the  orthogonal  projection  of  b  onto  R(A).  Thus,  the  LS 
vector  b'  can  be  viewed  as  the  perturbation  of  b  so  that  b'  can  be  generated  by  the  range  (or  the 
column  space)  of  A.  As  a  consequence,  the  LS  solution  assumes  that  errors  can  occur  only  in  the 
vector  b  and  not  in  the  matrix  A.  However,  this  assumption  is  not  always  realistic  as  errors  such 
as  sampling,  modeling,  and  instrumentation  errors  may  force  inaccuracies  in  A  also.  The  TLS  takes 
into  account  the  fact  that  errors  may  occur  both  in  the  matrix  A  and  in  the  vector  b  and  solves  the 
following  minimization  problem: 

minimize^  e  r«.(..d    l[A;h\  ~[A\h\  ll2 
subject  to    b_  e  R(A) 
Thus,  the  TLS  vector  6  can  be  viewed  as  resulting  from  the  perturbation  both  from  the  columns  of 


A  and  the  vector  b  so  that  6  belongs  to  R(A).  An  illustration  of  the  differences  between  the  LS  and 
the  TLS  solution  for  a  system  with  column  space  of  dimension  two  is  shown  in  Figure  1. 

22  SVD-Based  Solution  of  the  TLS  Problem 


The  Singular  Value  Decomposition  (SVD)  is  used  to  solve  the  TLS  problem  [1,2].    Recall 
that  the  system  to  be  solved  is  of  the  form  Ax  =  b,  which  can  be  reformulated  as: 


m 


-l 


=Q. 


The  TLS  solution  is  found  from  the  SVD  of  the  augmented  matrix 

C  =  [A;b]  =  U2VH. 


The  solution  is  said  to  be  generic  or  non  generic  depending  on  the  numerical  characteristics  of  the 
eigenvalue  matrix  E  and  the  right  singular  vector  matrix  V. 

2.2.a  TLS  Generic  Solution 


The  TLS  solution  xTLS  is  unique  and  said  to  be  generic  if  the  singular  value  on>an  +  l  and 
vn+in  +  i't0.   In  such  a  case  the  solution  is  given  by: 

X-TLS   =    •VVn+l.n  +  l[Vl,n  +  l VM+l]'- 

When  the  p  smallest  singular  values  are  equal,  i.e.,  o  >a  x  =  ...  =  on  +  1,  the  TLS  solution 
to  the  linear  system  is  not  unique,  and  any  linear  combination  of  the  singular  vectors  associated  with 
the  multiple  minimum  singular  value  can  be  chosen  provided  that  it  is  normalized  properly.  The 
resulting  minimum  norm  solution  can  be  shown  to  be  equal  to  [1,3] 


TLS 


V* 


C'c 


Matrices  V2  and  V2  are  submatrices  of  the  matrices  of  the  right  singular  vector  matrices  Vj  and  V2, 
where  V2  is  the  matrix  of  right  singular  vectors  associated  with  the  multiple  minimum  singular  value, 
and  Vj  contained  all  the  other  right  singular  vectors.  V1  and  V2  are  defined  as: 


y2  = 


and  Vx  = 


where  the  vectors  c1  and  g1  respectively  represent  the  first  row  of  V2  and  V2. 

Again,  the  TLS  solution  exists  only  if  all  vn+lj*0  for  i  =  n-p+  l,...,n+  1.   Van  Huffel  showed 
that  no  generic  solution  exists  only  when  the  matrix  A  is  nearly  rank  deficient,  or  when  the  set  of 


p+i 


:...=0 


n+1 


are 


equations  is  conflicting  [4].   The  set  of  equations  is  said  to  be  conflicting  when  on 

large.    In  such  a  case,  trying  to  model  the  data  using  a  linear  model  is  inaccurate,  and  a  better 


option  should  be  chosen  by  the  user.  When  <Vp+1=...=oll  +  1  are  small,  two  options  are  possible.  The 
first  option  for  the  user  is  to  remove  dependent  columns  of  A  to  reset  the  problem  so  that  the 
resulting  modified  A  matrix  is  regular,  and  to  solve  a  generic  TLS  problem.  How  to  pick  such 
columns  is  addressed  in  [2,14].  The  second  option  is  to  solve  a  nongeneric  TLS  problem  which  is 
obtained  from  adding  additional  constrains  on  the  generic  TLS  problem  in  order  for  it  to  be 
solvable. 

22.b  The  Nongeneric  TLS  Problem 

The  nongeneric  TLS  problem  is  addressed  in  details  in  Van  Huffel  et  al.  [4,5,7].  The 
resulting  solution  is  given  by: 

X-TLS    =    ■VVn  +  l.n-plVl.np'-"'Vn,n-pJ  ' 

when  on.p>  on.p+1  =  ....  on  +  1,  vn+1  .  =  0  for  i=  n-p+  1,...,  n+1,  and  vn  +  ln.p  *  0. 

23  Applications  of  the  TLS  to  Signal  Processing  Problems 

Numerous  researchers  in  Signal  Processing  have  reformulated  problems  in  terms  of  the  TLS 
technique.  Applications  can  be  found  in  array  processing  [11],  in  system  modeling  [10],  and  in 
frequency  estimation  of  sinusoidal  signals  [9,10-13]  for  example.  Results  show  that  improvements 
in  the  performances  can  be  obtained  when  reformulating  LS  problems  in  a  TLS  set-up  for  the 
applications  mentioned  above.  This  result  motivated  our  investigation  of  the  TLS  technique  to 
time-series  data  modeling.  ARMA  modeling  of  time-series  is  a  non-linear  problem  which  has  been 
extensively  studied  in  the  literature  [17].  The  Iterative  Prefiltering  (IP)  method  is  an  iterative 
linearized  formulation  which  was  originally  proposed  by  Steiglitz  and  McBride  to  identify  linear 
system  transfer  functions  [15,16],  and  the  IP  method  can  be  reformulated  to  model  time-series  data 
distorted  with  noise.  Note  that  the  method  is  not  insured  to  converge  to  the  optimal  solution  and 
De  Moor  showed  that  it  converges  to  a  suboptimal  solution  only  [18].  Nevertheless,  simulation 
results  have  shown  that  useful  models  can  still  be  obtained  using  the  IP  method  [19]. 

3.  The  Iterative  Prefiltering  Method 

3.1  Introduction 

The  Iterative  Prefiltering  method  attempts  to  model  a  time-series  data  with  an  ARMA(P,Q) 
system  using  a  time-domain  approach  by  minimizing  the  error  between  the  data  and  the  impulse 
response  of  the  system  to  be  estimated.   Using  the  Z-domain,  the  model  function  is  given  by: 

H(z)  -  *> 

A(Z) 


where 


p  p 

A(z)  =  J>(*)z"*    and    B&  =  £*(*)*"*• 

k=0  i=0 


The  error  function  for  the  problem  is  defined  by 


£(z)  =  xu-ito  -  wmzm. 

A{z)  A(z) 


(1) 


Non-linearity  in  the  problem  is  due  to  the  denominator  A(z)  in  the  error  function  E(z).  The 
problem  may  be  linearized  by  replacing  the  error  function  given  in  eq.  (1)  with  the  iterative  error 
function 


,(0,       .   X(zMw(z)-*w(z) 


E(,)(z)  = 


(i-i) 


C) 


where  indices  (i)  and  (i-1)  refer  to  iterations  (i)  and  (i-1).  At  each  iteration  the  quantities  A(l)(z) 
and  B(l)(z)  are  chosen  to  minimize  the  error  E(l)(z).  The  estimation  of  the  ARMA  transfer  function 
is  done  in  the  time  domain  using  the  time  series  x(n),  for  n  =  0,...,N-l,  by  rewriting  E(0(z)  in  the 
time-domain  as 


Eii)(n)=xli)H(n)*ali>-hli)H(n).biii.      n=0,...JJ-l 


(2) 


with: 

Eli>(n)=x{i)(n)M.a{i)-hm(ri).bii>    n=0,...^-l 

x}!!(n)  =  [x('\n\.,xi'\n-P)]T 

l£±(n)=[h('\n),.  ,h°\n-Q)]T 

a»=[l,a«(l),...,0W(P)] 

b^.=[b(i>(P),..,bl,\Q)] 

where  h(1)(n)  is  the  impulse  response  of  H0,(z)=  1/A(1)(z),  and  x(1)(n)  represents  the  output  of  x(n) 
through  the  filter  H(l)(z).    Eq.  (2)  expressed  in  matrix  form  becomes: 


£C0  = 


x(0(0) 
jc(0(0)        x(0(l) 


x('\N-P)    ...   x('\N-2)   x('\N-l) 


a(i)(P) 
a°\P-l) 

1 


/i(,)(0)         0 
h(,\l)      hu\0) 


h(,)(N-l) 


0 

b(,\0) 

0 

b(,\l) 

hm(N-Q) 

b('\Q) 

(3) 


which  may  be  rewritten  as: 


Equation  (4)  may  be  decomposed  into  two  parts  using  the  orthonormal  projection  matrix  P  onto 
the  column  space  of  X<0,  as  noted  by  McCleUan  and  Lee  [17],  which  leads  to: 

(5) 
with    P  =  H(i)[HmH{i)YlHm 

Recall  that  the  IP  iteration  minimizes  the  expression    ||E(,)||2,  which  leads  to  the  following 
minimization  problem  [17]: 

•  (a)  minimize  pxx(0o^  subject  to  a(i)(0)=  1, 

•  (b)  solve  ^=[HmH{i)YlHmX(i)(^. 

Step  (a)  above  may  be  rewritten  as: 
[C,-d][a(i»(P),...,l]T  =  0, 

where 

C=PXm[l:N,l:P],  d  =  PiX(i)[l:N,P+  1]. 

Step  (a)  of  the  IP  method  has  classically  solved  using  the  classical  LS  approach,  which 
accounts  for  errors  in  d  only.  However,  errors  occur  both  in  C  and  d.  A  better  fit  of  the  data  can 
be  obtained  by  taking  into  account  for  errors  in  both  in  C  and  d,  as  the  TLS  set-up  allows  us  to  do. 
Thus,  the  TLS-Based  IP  iteration  solves  for  (a)  using  a  TLS  approach  [20],  which  leads  to  the 
following  minimization  problem: 

•  (a)  Minimize  P^X(i)a(i)  subject  to  a(i)(0)=l,  using  a  TLS  method. 

•  (b)  solve  lP=[HmH{C)YxHmX^G^_. 
3.2  Data  Scaling 

3.2.a  Introduction 

Van  Huffel  has  shown  that  when  errors  in  the  TLS  matrix  and  right-hand-side  vector  are 
uncorrelated  with  zero-mean  and  equal  variance,  then  under  mild  conditions  the  TLS  solution  is 
a  strongly  consistent  estimate  of  the  true  solution  of  the  unperturbed  system  [4].  For  the  problem 
considered,  errors  in  the  data  result  from  errors  in  the  signal  x(n)  to  be  modeled.  Assuming  the 
data  to  be  modeled  is  distorted  with  Gaussian  white  noise,  the  errors  in  H(1)  are  correlated  as  they 
are  obtained  as  the  output  of  an  AR  linear  system.  However,  consistency  in  the  TLS  estimate  may 
still  be  retained  by  using  the  Generalized  TLS  (GTLS),  as  proposed  originally  by  Van  Huffel 
[2,4,17].  Thus,  data  scaling  is  needed  to  insure  diagonal  error  covariance  matrices,  and  the  problem 
becomes  to  estimate  the  error  covariance  matrices  F(l)  and  G(l)  at  each  iteration  given  by: 

F(i)  =  E[NH(i)N(i)], 
G(i)  =  E[N(i)N(i)H], 


where  N0)  represents  the  errors  (noise  contribution)  in  [C,-d]  at  iteration  (i). 

3.2.D  Computation  of  the  Matrix  F"1 

Assume  that  the  noisy  signal  x(n)  is  defined  by  x(n)  =  s(n)  +  w(n),  where  w(n)  is  wss  white 
noise,  and  s(n)  is  an  ARMA  signal.  The  noise  correlation  matrix  F("  may  be  expressed  in  terms 
of  the  correlation  matrix  Rw(l)  obtained  from  the  AR  process  1/A(1)(z): 

=  E[W®HPlHPi-W®'\ 

=  £[H^0///>xW/<0]  =  E[Wif)H(I-P)Wif)].  (6) 

=  E[  w^W®]  -E[  W®HW(i)] 

Computations  show  that  the  components  of  Ew(1)  in  equation  (6)  can  be  expressed  as: 

*)=EV^+«-^  (7) 

where  P,j  is  the  (i,j)lh  component  of  the  projection  matrix  P  defined  earlier,  and  rw(1)(n)  is  the 
correlation  function  obtained  when  passing  the  white  noise  w(n)  through  the  AR  system  with 
transfer  function  1/A(1)(z).  Using  the  fact  that  the  noise  distortion  w(n)  is  assumed  to  be  white, 
rw(l)(  |l  +  q-j-k| )  is  non  zero  only  when  l  +  q-j-k  =  0.  Thus,  the  double  summation  in  equation  (7) 
reduces  to  a  single  summation.  Furthermore,  the  correlation  function  rw(,)(k)  can  easily  be 
computed  using  the  results  by  Dugre  et.  al.  [6]  who  proposed  an  algorithm  to  generate  a  covariance 
sequence  from  its  AR  coefficients. 


3.2.C  Computation  of  the  matrix  G'" 

Again  the  error  covariance  matrix  G("  may  be  expressed  in  terms  of  the  correlation  matrix 
Rw(l)  obtained  from  the  AR  process  1/A(l)(z)  obtained  at  iteration  (i)  by: 

jpffl  =  E[N(f)N(i)H)  =  P±E[W{[)W{i)H)]P± 

(8) 

=  P'R§PX 


32.d  GTLS  Solution 

Applying  the  results  presented  by  Van  Huffel,  the  general  TLS-based  IP  iterative  solution 
for  a(l)  at  iteration  (i)  is  then  given  by  solving: 


(RTiC,-dKR?rl) 


kl\P) 


R1fO 


a(,\l) 


■Q    with    Rf=chol(F%    R$  =  U2*UH, 


with: 


H+=diag[Jox ,/^,0,...0] 


where  r  is  the  rank  of  RG(,).  Next,  the  estimate  for  b(l)  can  be  obtained  by  replacing  a(l)  by  its  value 
in  the  expression: 

I^_=[HmH^-lHmX^(^ 


33  Simulation  Results 

The  TLS-based  IP  method  is  implemented  on  time-series  data  generated  by  an  ARMA(3,4) 
in  additive  white  noise  and  its  performances  compared  with  that  of  the  LS-based  IP  method.  In  both 
cases,  initial  estimates  for  the  iterative  procedures  are  obtained  using  Yule  Walker  equations  for 
the  poles  and  Shank's  technique  when  solving  for  the  zeros.  The  resulting  noisy  signal  is  modelled 
with  an  ARMA(5,6)  using  a  sequence  of  length  60.  The  noise  variance  is  chosen  equal  to  2  and  the 
ARMA  parameters  are  chosen  equal  to: 

ap=  [1,-2,1.6725,-0.4613] 
bp=[l,-.9,0.04,0.3960,-0.2912] 

Figure  2. a  shows  the  noise-free  signal,  the  noisy  signal,  and  the  estimated  signal  obtained 
during  the  first  5  iterations  using  the  classic  IP  method.  Power  spectral  estimates  obtained  for  the 
noise-free  signal,  the  noisy  signal,  and  that  obtained  with  the  initial  ARMA  coefficients  estimated 
using  Yule  Walker  equations  and  Shank's  method,  are  shown  in  Figure  2.b.  Figures  3.a  and  3.b 
represent  the  same  estimates  obtained  when  using  the  proposed  TLS-based  IP  scheme.  Note  that 
Figures  2  and  3  represent  the  best  performances  obtained  during  the  first  5  iterations  for  both 
schemes.  Figure  4  and  5  represent  the  true  pole  (a  )  and  zero  (b  )  locations,  and  estimated  poles 
(a,  )  and  zeros  (b,  )  locations  obtained  when  using  the  IP  and  TLS-based  IP  schemes.  This  example 
illustrates  the  fact  that  the  bias  in  the  estimated  pole  locations  obtained  using  the  TLS-based 
technique  is  usually  smaller  than  that  obtained  using  the  LS-based  IP  method  when  the  SNR  is 
medium  to  high.  However,  we  noted  that  when  the  SNR  is  low,  no  distinct  improvement  is  noted 
consistently  when  using  the  proposed  technique. 

4.  Conclusions 


This  study  has  investigated  the  application  of  the  TLS  problem  in  the  Iterative  Prefiltering 
method  for  time-domain  ARMA  modeling.  Results  show  that  the  pole  bias  is  usually  smaller  when 
using  the  TLS-based  scheme  than  when  using  the  classical  LS-based  IP  method,  when  SNR  are 
medium  to  high.  We  note  that,  similarly  to  the  LS-based  method,  the  TLS-based  IP  method  is  not 


TLS-based  IP  method.     The  main  drawback  in  the  proposed  TLS-based  IP  method  is  the 
computational  load  increase. 
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(?ia; 


Figure  l.a)  LS  solution  to  Ax  =  b;  b)  TLS  solution  to  Ax  =  b. 
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Figure  2. a)  Noise-free  signal,  noisy  signal,  LS-based  IP  estimated  signal;  b)Spectra:  H  noise  free: 
noise-free  signal,  H_est: LS-based  IP  signal,  H_itl:  initial  estimate. 
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estimated  signal:  dash   -   noise  free  signal : solid  -  noisy  signal:  do  t 
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Figure  3. a)    Noise-free    signal,    noisy    signal,    TLS-based    IP    estimated    signal;    b)Spectra: 
H  noise_free:  noise-free  signal,  H_est:  TLS-based  IP  signal,  H_itl:  initial  estimate. 
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Figure  4.True  pole  (ap)  and  zero  (bp)  locations;  LS-based  IP  estimated  pole  (aip)  and  zero  (bip) 
locations. 
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Figure  5.True  pole  (ap)  and  zero  (bp)  locations;  TLS-based  IP  estimated  pole  (aip)  and  zero  (bip) 
locations. 
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